Interpretation of the expansion law of planetary nebulae 



Hiroko Matsumoto, Tsubasa Fukue, and Hideyuki Kamaya 
Department of Astronomy, School of Science, Kyoto University, Sakyo-Ku, Kyoto 606-8502 

kamaya@kusastro.kyoto-u.ac.jp 

(Received (reception date); accepted (acception date)) 
Abstract 

§ I We reproduce the expansion velocity-radius (Kxp^-Rn) relation in planetary neb- 

^ I ulae by considering a simple dynamical model, in order to investigate the dynamical 

evolution and formation of planetary nebulae. In our model, the planetary nebula is 
^ ■ formed and evolving by interaction of a fast wind from the central star with a slow 

wind from its progenitor, the AGB star. In particular, taking account of the mass loss 
history of the AGB star makes us succeed in the reproduction of the observed T4xp--Rn 
>■ ' sequence. As a result, examining the ensemble of the observational and theoretical 

^ I evolution models of PNe, we find that if the AGB star pulsates and its mass loss rate 

^ I changes with time (from ~ IO'^'^Mq jt^^ to ~ 10^^ Mq yr~^), the model agrees with 

O ! the observations. In terms of observation, we suggest that there are few planetary 

O ■ nebulae with larger expansion velocity and smaller radius because the evolutionary 

^ ■ time-scale of such nebulae is so short and the size of nebulae is so compact that it is 

I ' difficult for us to observe them. 

i-i • Key words: ISM: evolution — ISM: planetary nebulae: general — hydrodynamics 

> 

^ ' 1. Introduction 



Planetary nebulae (PNe) are the objects apart from the asymptotic giant branch (AGB), 
evolving to white dwarfs. It is thought that a stellar wind is important for the formation and 
evolution of PNe. To explain the features of PNe, Kwok et al. (1978) propose the interact- 
ing stellar wind (ISW) theory. After that many authors have studied PNe theoretically and 
observationally based on the ISW theory. 

Observationally, before the ISW theory was presented, a trend of the expansion velocity 
{Vcxp) of the gas of PNe was found (Kwok et al. 1978), and it has been interpreted in two 
different ways: one is in terms of a dynamical evolution of PN (Smith 1969; Bohuski and Smith 
1974; Sabbadin and Hamzaoglu 1982; Robinson et al. 1982), another is in terms of the evolution 
of a central star (Renzini 1979). 

Based on such studies, Sabbadin et al. (1984; hereinafter SGBO) proposed two models 
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to interpret that trend, the "two-wind" model and the "two-phase" model. The former means 
that the planetary nebula is the region of interaction of a fast wind emitted by the nucleus with 
a slow wind expelled by its progenitor. They added the effect of the momentum transfer from 
radiation field to the model suggested by Kwok et al.(1978) and Kwok(1982). The latter means 
that the planetary nebula is a Stromgren sphere evolving in an expanding nebula which results 
from a sudden ejection by the AGB star. They conclude only the "two-phase" model agrees 
with the observational data because the "two- wind" model can not reproduce the observed mass 
(Mn) and radius (Rn) relation. In the two-wind model the resulting mass of nebula is much 
smaller than that derived from observations (the mass decrement problem), while it matches 
the Ve^p-Rn relation. The mass decrement problem of the two-wind model is also stressed and 
examined in Schmidt- Voigt & Koppen (1987a,b). 

The "two-wind" model had another fault. SGBO didn't take account of the effect 
of shock sufficiently. As Marten & Schonberner (1991) stated, the simplification of SGBO 
corresponds to the fact that the width of the shell is infinitesimally thin. To resolve the mass 
decrement problem, it is essential to examine the width of the shell (Marten & Schonberner 
1991). Then, we confirm in this paper that it is easy to resolve the mass decrement problem if 
the width of the shell is determined as the effect of the shock. In this standing point, we try to 
study the two-wind model to satisfy M^-R^ and V^xp-Rn relations. We would fikc to comment 
that after SGBO, they developed their study on PNe furthermore (e.g. Sabbadin et al. 2005). 

We also examine the effect of the time-dependent mass loss rate of the AGB star, then we 
succeed in reproduction of the observed mass-radius relation simultaneously with the observed 
Vexp-Rn relation. The importance of the mass-loss to explain the evolution of PNe has been 
stressed by Marten & Schonberner (1991). In this paper, we insist furthermore that our modified 
"two- wind" model reproduces the ensemble of observational data adequately if the radial density 
profile of AGB matter has a sudden decrement at about O.lpc from the central star. This 
decrement is the result of the time-dependent mass-loss history of the AGB star. Thus, we 
reformulate the model concerning the mass-loss history of the AGB star. 

Thus, to explore the evolution and formation of PNe, we revisit the "two-wind" model 
in this paper. In §2 we describe our model. In §3 the results of the calculation compared with 
the observational data are presented. In §4, we discuss the results, and summarize our study 
in §5. 

2. Model description 

According to SGBO, there are two problems in the classical "two-wind" model: (a) the 
resulting mass of the PN is much smaller than that determined from observations [hereinafter 
problem(a)], (b) the "two- wind" model can't reproduce the observed Rn-M^ relation [See Fig.6 
and Fig. 7 of SGBO and hereinafter problem (b)]. Then, we reformulate the problem here. 
According to Frank (1994), it is important for us to study the evolution of PNe as a problem 
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Fig. 1. The evolutionary tracks of central stars. Horizontal axis is effective temperature, and vertical axis 
denotes luminosity. In this figure, three different lines correspond different mass: O.SSSM© (dotted line), 
O.6O5M0(sohd line), O.565M0 (dashed line). 



of radiation hydrodynamics. We admit it. However, if we are concerning the ensemble of the 
evolution of PNe, some simple models of PN evolution is useful, since the overall property 
of evolution of PNe is depicted as a combination among simple physics. Then, we try to 
reconstruct the two-wind model and explore its possibility as a tool to study PN evolution. 

2.1. General description 

Our model is based on the "two-wind" model which is proposed by SGBO. In the "two- 
wind'' model, the authors assume that PN is the region of interaction of a fast wind from the 
central star with a slow wind from the AGB star, considering the effect of a momentum transfer 
from the radiation field to the nebula motion. We assume a shell of PN is a shock region and 
the nebula evolution obeys the following equations: 

dM^ M^gb(0 ... 

~ -^n (1) 



and 



dt 

dVn 1 
dt ~ Mn 



« ^F^ K - Kgb ' + - - ^^n ' 



(2) 



where M^^it) and \4gb are the mass loss rate and the velocity of the wind from the AGB star; 
L, m, and v are the luminosity, the mass loss rate, and the velocity of the wind from the central 
star; Mn, i?n, and are the mass, the radius, and the expansion velocity of the PN; c is the 
light velocity; a is the efficiency factor of the transfered momentum. 

Eq.(l) means that the mass of PN is consisted of the total matter included in the radius 
i?n thanks to the shock. Magb(i) denotes the mass loss rate measured at the position of the 
expanding shell, and then it is represented of time explicitly. 
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Fig. 2. Mass-loss rate of a central star. Horizontal axis is evolutionary time [yr] and vertical axis is 
mass loss rate [Mq /yr] . Each of three lines corresponds to different central mass: dotted line corresponds 
0.836Mq, solid line is 0.605Mq, dashed line is O.SeSM©. 



In Eq.(2), the first term in the right hand side represents the contribution of the transfer 
of momentum from the radiation field to the nebula, the second term is the interaction with 
the AGB slow wind, and the third term corresponds to the interaction with the fast wind from 
the central star. 

2.2. Models of the central star 

Adopting Eqs.(l) and (2), wc take into account the evolution of luminosity L and the 
mass-loss m of the central star more realistically than SGBO adopted, while those evolution 
models are simplified as presented below. 

We assume that the mass of the central star is 0.605Mq because most nuclei of PNe are 
about 0.6Mq. Then we approximate the evolutionary track proposed by Blocker (1995b) shown 
in Fig.l. Furthermore, we compute the mass-loss rate shown in Fig. 2, according to Martin & 
Schonberner (1991). 

2.3. Initial Conditions 

As stated, SGBO insist there are two problems (mass decrement problem and R^-M^ 
relation problem) for the "two- wind" model. In order to solve these problems, we adopt Eq.(l) 
because the shock condition has to be satisfied there and the variation of the mass-loss rate 
of the AGB star to the "two-wind" model should be considered. Problem (a) is solved easily 
by only considering the shock condition. The previous model did not count the total mass 
reasonably. 

Following Blocker (1995a), we adopt the time-depending simplified mass-loss rate from 
the AGB star described in Fig. 3 to reproduce the i?n-Mn relation. In Fig. 3, the horizontal axis 
represents the time before the star enters post-AGB phase, i.e. the right edge of it is 1000 yr 
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Fig. 3. The mass-loss rate of the AGB star. Horizontal axis denotes time before a star enters PN phase, 
and vertical one is mass loss rate of the AGB star. 



before the post- AGB phase and the left edge corresponds 10^-^yr before it. For —t < 1000, Magb 
decreases slowly and connects rh in Fig. 2. 

On these bases, we set, as initial conditions, i?n=0.01pc, Mn=10~^-^Mo, i'n=Kgb+C's 
(Cgi the sound speed in pre-shock matter) at t—0, and Kigb=10 km s"^ which is kept during all 
the evolutionary time. These are chosen simply since PN emerges as a shocked layer. From t—0 
to ^=100000 yr, we calculate some models with different wind velocities from the central star: 
500 km s"^ < V < 10000 km s^^ at the intervals of 500 km s~^. The initial set of the parameters 
of Rn and are very consistent to those expected in a model calculation by Mellema (1994). 
According to the estimates of the surface brightness in Fig. 10 and the ionized mass in Fig. 19 
of Mellema (1994), when the PN shell is formed and fully ionized, the radius of the shell is less 
than 0.05 pc and its mass is less than 10~^ Mq. Thus, we confirm that we choose appropriate 
values as the initial conditions from the study of Mellema. 

2.4- Observational sample selection 

In our analysis, we use the [OIII] data listed in SGBO (See their table. 1 and Fig.l). 
We adopt this classical set of data because number of sample sets, in which both the size and 
expansion velocity are observed directly, does not increase even in recent years very much. 
Indeed, for our purpose, we must remove some objects whose expansion velocities are uncertain 
(small points in Fig.l of SGBO). 

According to Schonberner, Jacob, & Steffen (2005), Schonberner, et al. (2005), and so 
on, not only there are different sizes correspond to different components, optical imaging usually 
measures the size of the ionization front, not the actual size of the shell. Their statement is 
important if we are concerning the detailed structure of PNe. Fortunately, however, we are 
interested in the ensemble of the evolution of PNe, and do not examine the detailed structure 
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of PNe. Furthermore, the shell can be almost ionized after about 1000 years (e.g. Perinotto et 
al. 2004). Then, we expect that as long as the overall trend of PNe is concerned, our simple 
approach is still useful. 

3. Results 

3.1. Characteristic cases 

The results of our calculation in the case of Cs=5km s~^ and q;=0.5 are shown in Fig. 4, 
Fig. 5, and Fig.6. 




Fig. 4. Rn-Mn relation. Horizontal axis is radius in the unit of pc, vertical one is mass in the unit of 
solar mass. The bottom line is v = 500 km s~^, and the top line is = 10000 km s~^. 




Fig.4 represents Rn-M-^ relation of our calculation. From this figure, because of the 
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Fig. 6. 2Vexp-Rn relations. Horizontal axis denotes expansion velocity [km s~^], vertical one is radius of 
PN [pc]. The most left dashed line corresponds w=500 [km s~^], the most right one is w=10000 [km s"-'^]. 
Plus sign represents observational data, cross corresponds the point which evolutionary time is ^=1000 yr, 
and asterisks are marked from t=2000 yr at the intervals of 2000 yr. 



time-depending AGB mass-loss rate, the increasing rate of Mn decreases after PN has size of 
i?n ~0.1pc, then the increase rate turns to be smalL Thus, we can succeed in the reproduction 
of the observed tendency, and resolve the problem (b). 

Around the right top of Fig.4, Mn increases rapidly again. This is because the pulsation 
of the AGB star finishes (See Fig. 3, it corresponds log{—t) ~ 4.5). However it doesn't affect the 
value of expansion velocity up to i?n=0.3 pc. It influences only the resulting mass of PN. Fig.5 
is time-mass plane. The bottom line is the model i;=500 km s^^, the top is the model f =10000 
km s^^. From our simple models, it takes over 6000 yr that the radius of model nebula reaches 
0.3 pc (See Fig.6), and the resulting mass is over O.IMq when i?n ~ 0.3 pc, thus we can solve 
the problem (a) in the same time. 

In Fig.6, the relation of 2Vcxp-Iiu is presented as one of main results in our study. We 
can almost re-product the distribution of data points by changing only the fast wind speed. 
The most left curve is the model v =500 km s~^ and the most right one corresponds v =10000 
km s"-*^. Prom Eq.(2), if v has large value, the shell of PN is accelerated strongly in the early 
phase then has a large expansion velocity. The expansion velocity is decelerated as the PN 
evolutes, because of the decrease of the luminosity and mass loss rate of the central star. We 
also over-plot the evolutionary time as asterisks from t=2000 yr at steps of 2000 yr. One can 
comprehend from it that most of data points exist in the region of t=2000 yr-4000 yr. Star 
which mass is O.6O5M0 enters the realm of white dwarf at i ~ 7000 yr (See Fig.l), so PN phase 
must be before 7000 yr. Our results are consistent in this meaning. 
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Fig. 7. 2Vexp-Rn:Cs=5-10 km s ^. Horizontal axis is expansion velocity [km s and vertical one is 
radius [pc]. Plus sign represents observational data. 



3.2. Parameter dependence 

Setting the fast wind velocity f=1000 km s~^, we inquire into the effects of other pa- 
rameters: the sound speed (Cs), the transfer efficiency (a), and the slow wind velocity (V^agb)- 
The effect of other parameters to the results is essentially the same as examined in SGBO. 

3.2.1. Sound speed 

Fig.7 shows the results of different sound speed in the pre-shock matter: Cs=5-10 km 
s~^ at the intervals of 0.5km s~^. Even if Cg changes, i.e. the temperature in the pre-shock 
matter changes, V^^-R^ relation is not so different. The variation of Cg affects only Rn-M^ 
relation, as seen in Fig. 8. It is because we set that the mass of PN increases along Eq.(l) only 
as long as the shock occurs. The larger the sound speed is, the shorter the period during which 
the shock is occurring is and the smaller the final mass of PN is. If the shock doesn't occur, 
the mass of PN has a constant value in our model. 

3.2.2. Transfer efficiency 

The sequences of different transfer efficiency are shown in Fig. 9 and Fig. 10. We vary 
q;=0-1 at the steps of 0.1. As a increases, the influence on the evolution of PN due to the 
variation of central star's luminosity becomes large. 

3.2.3. Slow wind velocity 

The 2V^xp--Rn relations are drawn in Fig. 11, and the R^-M^ relations are in Fig. 12 when 
we change the slow wind velocity from 5km s^^ to 24km s^^ at the intervals of 1km s~^. In 
Fig. 11, the most left line corresponds K,gb=5km s^^ and the most right one does Vjjgb=24km 
s"^. Also in Fig. 12, the top curve is in the case of \4gb=5km s~^ and the bottom is Kigb=24km 
s~^. If we allow a turning point, at which the increasing rate of mass turns to be small, to have 
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Fig. 8. log Rn-log Mn:Cs=5-10 km s ^. Horizontal axis corresponds radius, and vertical axis is mass. 
The bottom line is in the case of Cs=10 km s~^. 
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Fig. 9. 2Vexp-Rn:alpha=0-l. Horizontal axis is expansion velocity [km s~^], and vertical axis denotes 
radius [pc]. The most left dashed line is a=0, the most right one is a=l, and plus sign means observational 
data. 



uncertainty of a factor 3 in log (i?n)-log (Mn) plane, ones can find that the difference of AGB 
wind velocity also explains the distribution of observed data well. However, AGB wind velocity 
is not so large (<30 km s~^) that the models varying only \4gb does not cover the region with 
large expansion velocity. Prom this point of view, we conclude the central star's wind velocity 
is more proper parameter which explains 2\4xp--Rn relation. 

4. Discussion 

In this section, we examine our model compared with the basic model basing on SGBO 
and discuss the physics of PN evolution. 
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Fig. 10. log Rn-log Mn:alpha=0-1. Horizontal axis is radius, and vertical axis is mass. The bottom line 
is in the case of Cs=0. 




Fig. 11. 2Vexp-Rn:Vagb=5-24km s~^. Horizontal axis denotes expansion velocity [km s~^] and vertical 
axis is radius [pc]. Plus sign corresponds observed data, the most left dashed line is in the case of V^gb=5km 
s~^, and the most right one is the model with ^gb=24km s~^. 



4-1. Difference from a basic model 

As already mentioned, there are two problems in the classical two- wind model by SGBO: 
(a) the decrement of resulting mass, and (b) the discrepancy in the R^^-Mj^ relation between 
observations and the model. Although SGBO's model explains the 2\4xp--Rn relation, the 
authors abandon their idea because of the problems (a) and (b). On the other hand, we take 
account of the effect of shock and time-depending mass loss rate of the AGB star, then the 
problems (a) and (b) are solved easily. Prom this point, we convince for the evolution of PN to 
consider the influence of shock and variation of mass loss rate of the AGB star. That is, the 
AGB star is pulsating so its mass loss rate varies with time during AGB phase, and during the 
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PN phase the fast wind speed is much larger than the slow wind speed to be a shock: v » T4gb- 
4-2. Physical origin of the Vexp-Rn relation 
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Fig. 13. 2Vexp-Rn:Magb=2xlO~^M0 yr~^. Horizontal axis is expansion velocity [km s~^], and vertical 
axis is radius [pc]. Plus sign means the data. Magb=2xlO~^M0 yr~^, a=0.5, Cs=10 km s~^ and other 
parameter's value are written in the figure. 



Based on Sect. 4.1, we discuss physical background of PN evolution in this subsection. In 
our model to explain the 2Ve^p-R^ relation, the greatest point is considering the time-depending 
mass loss rate of the AGB star. We also calculate the models which Magb value is a constant. 
Some results are presented in Fig. 13, Fig. 14, and Fig. 15. The sequences of Magb=2xlO~^M0 
yr"^ are in Fig.13, those of Magb=3xlO~^MQ yr"^ are in Fig.l4, and those of Magb=4xlO-^M0 
yr~^ are in Fig. 15. In these calculations, we set a— 0.5, Cs=10 km s~^ and examine the range 
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Fig. 14. 2Vexp-Rn:Magb=3xlO-5MQ yr'^ Same as Fig.l3 but in this figure Mlgb=3x IQ-^Mq yr'^ 
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Fig. 15. 2Vexp-Rn:Magb=4xlO-5M0 yr'^ Same as Fig.13 but in this figure Mlgb=4xlO-5MQ yr"!. 



of ^=500-10000 km s^^ at the steps of 500 km s^^ The three cases of V;gb=10,15,20 km s"^ 
are examined. In those figures, we plot only the models which satisfy the R^-M^ relation if the 
lines turn around the point of log(Rn) ~ -1.0 and log(Mn) ~ -1.0 in log(Rn)-log(Mn) plane. The 
R-o-M^ relations of the selected models are shown in Fig. 16, Fig. 17, and Fig. 18. 

From Figs. 13, 14, and 15, when Magb=const, we can re-product only the region of 2V^xp < 
70 km s~^. The models which are able to trace the realm of larger expansion velocity do not 
match the R^^-M^ data. 

Comparing time-depending mass loss models (Fig. 6) with steady mass loss models 
(Figs. 13, 14, 15), one can know how important for the evolution of PNe to take account of 
the variation of AGB mass-loss rate. In other words, the radial density distribution of ambient 
matter around PN, p(r), is not in proportion to r^^. It needs to fall suddenly at r ~0.1pc due to 
the pulsation history during AGB phase. Furthermore, this indicates the difference between the 
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Fig. 16. log Rn-log Mn:Magb=2x 1O~^M0 yr~^. Horizontal axis denotes radius and vertical axis is mass. 
-^agb=2xlO~^M0 yr~^, o;=0.5, Cs=10 km s~^ and other parameter's value are written in the figure. 
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Fig. 17. log Rn-log Mn:Mlgb=3xlO-^M0 yr'^ Same as Fig.l6 but Magb=3xlO-^M0 yr'^ 



dynamical age and evolutionary age of PNc since the expansion law is affected by the density 
profile around the proto-PNe. We expect the age-discrepancy problem of PNe (Mellema 1994) 
may be partially resolved if the origin of the velocity-radius relation is resolved. 

By the way, the variety of v contributes greatly to the explanation of data distribution 
in 2T4xp--Rn plane. In some models, v has a very large value up to 10000 km s~^ at t = 0. 
However, v may not attain such a large speed in the theory of line-driven wind. We deals with 
this problem in §§4.4 again. 

4-3. Observational implications 

In 2Vexp-Rn plane, most of observational data exist in the range of 2V^xp=40-80 km 
s~^ and it!n=0.05-0.15pc. We present by our calculations that these data are in the realm of 
t=2000-4000 yr(cf. Fig. 6 in Sect. 2) and it coincides with the fact that FN is an object between 

13 



the AGB star and white dwarf. Our results suggest the reason why only few data is in the 
smaller radius region. PNe evolve very quickly across this region, and then it is difficult for us 
to observe. In the larger radius region, it is also difficult to observe since PNe expand too large 
and diffuse. By a more highly accurate survey of PNe, in future, we will find that there are a 
lot of faint PNe and some of them have larger radius among them. 

4-4- Velocity problem 

As stated in subsection 4.2, we have a time gap problem concerning the fast wind velocity 
from the central star of PN. Along the line-driven wind theory, it takes about 8000 yr until v 
reaches 10000 km s~^, while in our calculation it needs to be attained at t=0. However we will 
solve this problem in two different simple ways: (1) accretion jet and (2) compactness of the 
central star. We shall examine the case of (1). In the hne-driven wind it is assumed the PN field 
is spherical symmetry. The assumption is not so realistic and many PNe have an asymmetric 
jet due to its magnetic field. In this way, the gas of PN is accelerated so quickly that our 
supposition for v is good. By this way of thinking, we may say that the reason why few PNe 
are in the larger expansion velocity and smaller radius region in their early phase (<2000 yr) 
is that the total number of asymmetrically accelerated PNe is small in their early evolution. 

The other explanation of the case of (2) is rather simple, although the acceleration 
mechanism can not be specified. If the central star has very small radius (~1000 km), its 
gravity at the surface of the star is so large. To escape from the central star gravity field, a 
matter must have a very large velocity. In fact, a radius of typical white dwarf is in the range 
of 1000-10000 km. It is expected that the acceleration mechanism of PNe will be proved by 
future observations. 

5. Summary 

Constructing the modified "two-wind" model, we find that the observational data are 
reproduced adequately if the radial density of the AGB matter has a sudden fall at about O.lpc 
from the central star. It corresponds the AGB star pulsates and its mass loss rate increases. 
When we set the start point of evolutionary track of PN is t—0 and the velocity of the AGB 
matter (Kgb) is 10 km s'^, we estimate Magb=lO-^'^M0 at f=-10^ Vr to M^g^=10-^MQ at 
t——10^'^jr. It means that the features of each PN considerably depend on its progenitor's 
mass loss history, and by only including the effect of time-dependent mass loss rate of the AGB 
star we can explain the observed data in the context of fluid dynamics. Also our model gives 
an observational suggestion that the PNe with larger Vexp (> 25km s~^) and smaller radius(< 
O.lpc) are evolving so quickly that it is difficult to observe them. However, to completely 
account the observational trend, we should investigate the acceleration mechanism of the fast 
wind since as long as spherical symmetry model is adopted, we need the long time-scale to 
reach at 10000 km s~^. In order to prove the acceleration mechanism of the shell of PN, we 
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have to do a more accurate survey of PNe because observations of only very young PNe make 
us know the mechanism. We expect that such observations will be done in the near future. 
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